library(cregg)
library(ggplot2)
library(texreg)

d <- read.csv("CY_conjoint_replication_data.csv", stringsAsFactors=T)
names(d)

#### data management ####
# shorten labels 
d$territ <- as.character(d$terr)
d$territ[d$territ=='plus Morphou'] <- 'plus Morphou'
d$territ[d$territ=='but Morphou to stay in TC administration'] <- 'Morphou to stay in TC admin'
d$territ[d$territ=='Morphou and North Karpasia to become federal areas'] <- 'Morphou and N.Karpasia federal areas'
d$territ[d$territ=='plus Morphou, Karpasia, Yialousa'] <- 'plus Morphou, Karpasia, Yialousa'
d$territ[d$territ=='plus old part of Morphou, Karpasia, Yialousa'] <- 'plus old part of Morphou, Karpasia, Yialousa'
d$territ <- as.factor(d$territ)

d$exec <- as.character(d$fed)
d$exec[d$exec=='support of at least a quarter of MPs from each community'] <- 'support of at least 1/4 of MPs'
d$exec[d$exec=='a majority in the assembly or voters regardless of ethnicity'] <- 'a majority in the assembly'
d$exec[d$exec=='all parties in proportion to their seats in the assembly'] <- 'all parties in proportion to their seats'
d$exec[d$exec=='GC and TC co-presidents elected through cross-voting'] <- 'GC and TC co-presidents'
d$exec[d$exec=='GC president and TC vice-president elected through cross-voting'] <- 'GC president and TC vice-president'
d$exec <- as.factor(d$exec)

d$implem <- as.character(d$impl)
d$implem[d$implem=='UN with the three former guarantors Greece, Turkey and the United Kingdom'] <- 'UN with the three former guarantors'
d$implem[d$implem=='UN with EU countries such as Ireland, France and Germany'] <- 'UN with EU countries'
d$implem[d$implem=='UN with third countries such as Japan, Australia and Canada'] <- 'UN with third countries'
d$implem[d$implem=='UN with third party such as NATO'] <- 'UN with third party'
d$implem <- as.factor(d$implem)

d$courts <- as.character(d$sup)
d$courts[d$courts=='with equal number of GCs & TCs with a minority of judges appointed by the ECHR'] <- 'equal no. of GCs & TCs'
d$courts[d$sup=='by a special international UN tribunal with headquarters in Cyprus'] <- 'by a special international UN tribunal'
d$courts[d$sup=='with equal number of GCs & TCs with rotating chair'] <- 'equal no. with rotating chair'
d$courts[d$sup=='with a majority of judges appointed by the ECHR'] <- 'a judges majority appointed by the ECHR'
d$courts <- as.factor(d$courts)
table(d$courts)

# relevel
levels(d$fed)
d$exec <- relevel(d$exec, ref =  "support of at least 1/4 of MPs")
d$exec <- droplevels(d$exec)
d$territ <- relevel(d$territ, ref = "plus Morphou")
d$proper <- relevel(d$proper, ref =  "50,000 Euro")
d$courts <- relevel(d$courts, ref =  "equal no. of GCs & TCs")
d$implem <- relevel(d$implem, ref = "UN with the three former guarantors")

# set feature labels
attr(d$exec, "label") <- "EXECUTIVE"
attr(d$territ, "label") <- "TERRITORIAL ISSUES"
attr(d$proper, "label") <- "ECONOMIC RESTITUTION"
attr(d$courts, "label") <- "COURTS"
attr(d$implem, "label") <- "IMPLEMENTATION" # should this be called legislative?
attr(d$agegroup, "label") <- "AGE"
attr(d$education, "label") <- "EDUCATION"
attr(d$gender, "label") <- "SEX"
attr(d$Income, "label") <- "INCOME"

#### AMCE ####
# without weights
# formula model conjoint
m <- selected ~ exec + territ + proper + courts + implem  

amces <- cj(d, m, id = ~IDM,  by = ~ community, feature_order = c('territ', 'exec',  'proper', 'courts', 'implem'))

head(amces[c("feature", "level", "estimate", "std.error")], 10L)  
# plot 
# black for TC and orange for GC
plot(amces, group = 'community',  header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 


#### MARGINAL MEANS ####
marginalmeans <- cj(d, m, id = ~IDM,  by = ~ community, estimate = 'mm')
head(marginalmeans[c("feature", "level", "estimate", "std.error")], 20L) 


# plot 
plot(marginalmeans, group = 'community', vline = 0.5,  header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

# APPENDIX
#### Regression table for AMCE ####
# GC
m_amce_gc <- lm(selected ~ territ + exec + courts + proper + implem, subset(d, community=='GC'))
summary(m_amce_gc)

# TC
m_amce_tc <- lm(selected ~ territ + exec + courts + proper + implem  , subset(d, community=='TC'))
summary(m_amce_tc)

# create table 
wordreg(list(m_amce_gc, m_amce_tc),
        file = "CY_amce.doc", 
        ci.force = TRUE,
        groups = list("TERRITORIAL ISSUES" = 2:5, "EXECUTIVE" = 6:9, "ECONOMIC RESTITUTION" = 10:13, "COURTS" = 14:16, "IMPLEMENTATION" = 17:19),
        custom.coef.names = c("(Intercept)", 
                              "but Morphou to stay in TC administration", "Morphou and North Karpasia to become federal areas", "plus Morphou, Karpasia, Yialousa", "plus old part of Morphou, Karpasia, Yialousa",
                              "a majority in the assembly or voters regardless of ethnicity", "all parties in proportion to their seats in the assembly", "GC and TC co-presidents elected through cross-voting", "GC president and TC vice-president elected through cross-voting",
                              "by a special international UN tribunal with headquarters in Cyprus", "with a majority of judges appointed by the ECHR", "with equal number of GCs & TCs with rotating chair",
                              "150,000 Euro", "200,000 Euro", "300,000 Euro", "300,000 Euro plus housing",
                              "UN with EU countries such as Ireland, France and Germany", "UN with third countries such as Japan, Australia and Canada", "UN with third party such as NATO"), 
        custom.model.names = c("GC","TC"))
unlink("CY_amce.doc")

#### Regression table for MM ####
mm <- marginalmeans[c("feature", "level", "estimate", "p", "lower", "upper", "community")]
# scientific is an issue!
write.table(format(mm, digits=2), file = "CY_mm.txt", sep = ";", quote = FALSE, row.names = F)

#### AMCE with weights ####
# with weights, by community 3 ways
amces_weights <- cj(d, m, id = ~IDM, weights = ~weights, by = ~ community, feature_order = c('territ', 'exec',  'proper', 'courts', 'implem'))

head(amces_weights[c("feature", "level", "estimate", "std.error")], 10L) 

# plot 
plot(amces_weights, group = 'community',  header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### MM with weights ####
marginalmeans_weights <- cj(d, m, id = ~IDM,  weights = ~weights, by = ~ community, estimate = 'mm')
head(marginalmeans_weights[c("feature", "level", "estimate", "std.error")], 20L) 

# SET THE SAME LIMITS TO THE X SCALE GO WITH 4 OR WITH 2?
# plot 
plot(marginalmeans_weights, group = 'community', vline = 0.5,  header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### AMCE with controls ####
# without weights
# drop the missing on age, edu and gender and unused levels
d <- d[!(d$agegroup==""), ]
d$agegroup <- droplevels(d$agegroup)
levels(d$agegroup)

table(d$education)
d <- d[!(d$education==""), ]
d$education <- droplevels(d$education)
levels(d$education)

d <- d[!(d$gender==""), ]
d$gender <- droplevels(d$gender)
levels(d$gender)

d <- d[!(d$Income==""), ]

d$Income <- droplevels(d$Income)
d$Income <- as.character(d$Income)
d$Income[d$Income=='High Inc'] <- 'High Income'
d$Income[d$Income=='Low Inco'] <- 'Low Income'
d$Income[d$Income=='Middle I'] <- 'Middle Income'

d$Income <- factor(d$Income, levels = c('Low Income', 'Middle Income', 'High Income'))
levels(d$Income)

# formula model conjoint with controls

m_controls <- selected ~  exec + territ + proper + courts + implem + agegroup + gender + education + Income
amces_controls <- cj(d, m_controls, id = ~IDM,  by = ~ community)

# plot 
plot(amces_controls, group = 'community',  header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 

#### Regression table for AMCE with controls without weights ####
# GC
m_amce_contr_gc <- lm(selected ~ territ + exec + courts + proper + implem + agegroup + gender + education, subset(d, community=='GC'))
summary(m_amce_contr_gc)

# TC
m_amce_contr_tc <- lm(selected ~ territ + exec + courts + proper + implem + agegroup + gender + education, subset(d, community=='TC'))
summary(m_amce_contr_tc)

# create table 
wordreg(list(m_amce_contr_gc, m_amce_contr_tc), 
        file = "CY_amce_controls.doc",
        ci.force = TRUE, 
        groups = list("TERRITORIAL ISSUES" = 2:5, "EXECUTIVE" = 6:9, "ECONOMIC RESTITUTION" = 10:13, "COURTS" = 14:16, "IMPLEMENTATION" = 17:19, "AGE" = 20:22, "GENDER" = 23, "EDUCATION" = 24:27),
        custom.coef.names = c("(Intercept)", 
                              "but Morphou to stay in TC administration", "Morphou and North Karpasia to become federal areas", "plus Morphou, Karpasia, Yialousa", "plus old part of Morphou, Karpasia, Yialousa", 
                              "a majority in the assembly or voters regardless of ethnicity", "all parties in proportion to their seats in the assembly", "GC and TC co-presidents elected through cross-voting", "GC president and TC vice-president elected through cross-voting",
                              "150,000 Euro", "200,000 Euro", "300,000 Euro", "300,000 Euro plus housing", 
                              "UN with EU countries such as Ireland, France and Germany", "UN with third countries such as Japan, Australia and Canada", "UN with third party such as NATO", 
                              "by a special international UN tribunal with headquarters in Cyprus", "with a majority of judges appointed by the ECHR", "with equal number of GCs & TCs with rotating chair", 
                              "18-34", "35-54",  "55+", "Male", "Elementary", "High School", "Postgraduate", "Secondary", ''),
        custom.model.names = c("GC","TC"))
unlink("CY_amce_controls.doc")

#### AMCE first conjoint only ####
d_first <- d[d$task==1,]
# set feature labels
attr(d_first$exec, "label") <- "EXECUTIVE"
attr(d_first$territ, "label") <- "TERRITORIAL ISSUES"
attr(d_first$proper, "label") <- "ECONOMIC RESTITUTION"
attr(d_first$courts, "label") <- "COURTS"
attr(d_first$implem, "label") <- "IMPLEMENTATION"

amces_first <- cj(d_first, m, id = ~IDM,  by = ~ community, feature_order = c('territ', 'exec',  'proper', 'courts', 'implem'))

plot(amces_first, group = 'community',  header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 


#### AMCE last conjoint only ####
d_last <- d[d$task==5,]
# set feature labels
attr(d_last$exec, "label") <- "EXECUTIVE"
attr(d_last$territ, "label") <- "TERRITORIAL ISSUES"
attr(d_last$proper, "label") <- "ECONOMIC RESTITUTION"
attr(d_last$courts, "label") <- "COURTS"
attr(d_last$implem, "label") <- "IMPLEMENTATION"

amces_last <- cj(d_last, m, id = ~IDM,  by = ~ community, feature_order = c('territ', 'exec',  'proper', 'courts', 'implem'))

plot(amces_last, group = 'community', header_fmt = "%s") +
  scale_color_manual(values = c("#FF6600", "#000000"), na.translate = F) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  theme(legend.title = element_blank(),
        strip.text.y = element_blank()) 


#### Omnibus F test ####
# cj_anova does not support weights
cj_diff <- cj_anova(d, m, id = ~IDM,  by = ~ community)
cj_diff 
# present in table
write.table(format(cj_diff), file = "CY_omnibus_f_test.txt", sep = ",", quote = FALSE, row.names = F)







